In this script we will look at some interactive plots showing the effect of sex on the face. Because this is a preliminary assessment, the effect of sex is estimated using the whole sample (without taking their ancestries into account). This script is built from the hypothesis testing in which I created the vectors of sexual dimorphism. Remember to open the file in chrome.

Preliminaries

We will load the necessary libraries, and datasets

library(ggplot2)
library(dplyr)
library(plotly)
library(tidyr)
setwd('..')
path <- getwd()
setwd(paste(path, "/Results/FacePCA", sep = ""))
coeffs   <- read.csv("coeffs.csv")
eigenvec <- read.csv("eigenvectors.csv", header=F)
eigenval <- read.csv("eigenvalues.csv", header=F)
means    <- read.csv("means.csv", header=F)
facets   <- read.csv("facets.csv", header=F)
setwd(paste(path, "/Results/FSD", sep = ""))
load("fitmodel.RData")
vectors <- read.csv("vectors.csv", row.names = 1)

Data cleaning

For data cleaning, we will remove every individual that is neither Male or Female, and individuals with NAs in any variable (Height, Weight, and Age).

coeffs = filter(coeffs, Sex == 'Female' | Sex == 'Male')
coeffs = droplevels(coeffs)
coeffs = na.omit(coeffs)
dim(coeffs)
## [1] 5374   93
head(coeffs[,1:8])
##      ID    Sex Age  Height Weight      BMI       PC1         PC2
## 1 12001 Female  21 164.200   74.4 27.59476 -2.749977 -0.04763979
## 2 12002   Male  19 162.400   78.6 29.80235 -3.001070  0.10574822
## 3 12004 Female  19 163.500   51.2 19.15290 -0.939990 -1.44340059
## 4 12005   Male  21 160.625   62.4 24.18568 -2.595005  1.39772255
## 5 12006   Male  21 182.500   76.5 22.96866  2.580508 -0.08551034
## 6 12007 Female  22 154.500   48.7 20.40196 -1.539039 -1.32810730

After cleaning, we end up with 5374 individuals

Face plots

In this section we will see the interactive plots showing facial sexual dimorphism. To do that we need to transform the PCA vectors to the shape coordinates. Each vector (total, allometric and non-allometric) is multiplied by a constant to exaggerate the facial effects (in this case 2 and -2), multiplied by the matrix of eigenvectors, and then added the means for each landmark.

fsdcoords <- data.frame(Tot.male = rep(0, nrow(means)), Tot.female = rep(0, nrow(means)), 
                        Allo.male = rep(0, nrow(means)), Allo.female = rep(0, nrow(means)),
                        Nallo.male = rep(0, nrow(means)), Nallo.female = rep(0, nrow(means)))
fsdcoords[,c(1,3,5)] <- apply(as.matrix(vectors[3:5,] * -2) %*% t(as.matrix(eigenvec)), 1, function(x) x + t(means))
fsdcoords[,c(2,4,6)] <- apply(as.matrix(vectors[3:5,] * 2) %*% t(as.matrix(eigenvec)), 1, function(x) x + t(means))
averagecoords <- t(( colMeans(vectors[1:2,]) %*% t(as.matrix(eigenvec))) + t(means))

The first set of plots show the overlaying of the two face shapes. In this case male faces are shown as turquoise, and females as wheat color.

source("PlotFaces.R")
Plot2Faces(fsdcoords[ ,1], fsdcoords[,2], facets, "Total SD")
Plot2Faces(fsdcoords[ ,3], fsdcoords[,4], facets, "Allometric SD")
Plot2Faces(fsdcoords[ ,5], fsdcoords[,6], facets, "Non allometric SD")
source("Distances.R")
eucdist <- data.frame(Totdist = rep(0, nrow(means)), 
                      Allodist =  rep(0, nrow(means)), 
                      Nallodist = rep(0, nrow(means)))

eucdist[,1] <- getDistance(fsdcoords[ ,1], fsdcoords[,2])
eucdist[,2] <- getDistance(fsdcoords[ ,3], fsdcoords[,4])
eucdist[,3] <- getDistance(fsdcoords[ ,5], fsdcoords[,6])

#On the same scale
eucdist <- gather(eucdist, type, dist)
eucdist$dist <- scales::rescale(eucdist$dist)

In this other set of plots you can see the euclidean distance between the two morphs, as a colormap in the average face. At this point, the colormap is not scaled to be comparable across plost, and the legend does not show correctly the euclidean distances.

source("PlotFaces.R")
Plot1Face(averagecoords, facets, colormap = c(as.matrix(eucdist %>% filter(type == "Totdist") %>% 
                                                          dplyr::select(dist))), title = "Total SD")
Plot1Face(averagecoords, facets, colormap = c(as.matrix(eucdist %>% filter(type == "Allodist") %>% 
                                                          dplyr::select(dist))), title = "Allometric SD")
Plot1Face(averagecoords, facets, colormap = c(as.matrix(eucdist %>% filter(type == "Nallodist") %>% 
                                                          dplyr::select(dist))), title = "Non allometric SD")